#include "error.inc"

MODULE GET_STL_DATA_MOD

   USE bc, only: bc_defined, bc_type, bc_type_enum, is_cg
   USE compar, only: mype, pe_io
   USE constant, only: pi
   USE cutcell, only: ijk_of_node, f_at
   USE cutcell, only: use_stl, use_msh
   USE error_manager
   USE is
   USE make_upper_case_mod, only: make_upper_case
   USE param, only: dimension_3
   USE param1, only: one, undefined, nl
   USE quadric, only: cross_product, n_quadric
   USE run, only: run_name
   USE stl
   use param1, only: zero, one
   use geometry, only: x_min,x_max,y_min,y_max,z_min,z_max,no_k
   use bc, only: bc_x_w, bc_x_e, bc_y_s, bc_y_n, bc_z_b, bc_z_t
   use DEFINE_QUADRICS_MOD

CONTAINS

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: GET_MSH_DATA                                           C
!  Purpose: reads face vertices and normal vectors from an MSH file    C
!           (generated by Gambit). After reading the geometry, the     C
!           data is converted into same format as STL and the          C
!           pre-processing switches to STL procedure                   C
!                                                                      C
!  Author: Jeff Dietiker                              Date: 30-JAN-09  C
!  Reviewer:                                          Date: **-***-**  C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

   SUBROUTINE GET_MSH_DATA

      IMPLICIT NONE

      INTEGER,PARAMETER :: MAX_POINTS = 10000000    ! Currently limited to 10 Million, increase if needed
      INTEGER,PARAMETER :: MAX_ZONES  = 1000        ! Currently limited to 1,000     , increase if needed

      INTEGER ::I,D,NN,NF

      INTEGER   :: GRID_DIMENSION, N_POINTS,N_FACES,N_FACE_ZONES
      INTEGER   :: ZONE,ZONE_ID,N1,N2
      INTEGER   :: PPFACE

      LOGICAL :: ALL_POINTS_READ,ALL_FACES_READ

      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::POINT_ZONE_INFO,FACE_ZONE_INFO
      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: POINT_COORDS

      CHARACTER(LEN=18), ALLOCATABLE, DIMENSION(:,:) :: BC_LABEL_TEXT

      LOGICAL, ALLOCATABLE, DIMENSION(:) :: BC_ASSIGNED

      CHARACTER(LEN=18),DIMENSION(10) :: BUFF_CHAR

      INTEGER :: P1,P2,P3,P4
      DOUBLE PRECISION, DIMENSION(3) :: NORMAL, VECTMP, VECTMP2
      DOUBLE PRECISION ::NORM

      DOUBLE PRECISION ::v1x,v1y,v1z
      DOUBLE PRECISION ::v2x,v2y,v2z
      DOUBLE PRECISION ::v3x,v3y,v3z
      DOUBLE PRECISION ::ABSTRANS

      LOGICAL :: PRESENT

      INTEGER :: BC_LABELS_TO_READ, BC_LABELS_READ,NFZ,ZID,L2
      INTEGER :: N_BC_IGNORED
      INTEGER :: N_QUAD2TRI
!---------------------------------------------------------------------

      ALLOCATE(POINT_COORDS(3,MAX_POINTS))
      ALLOCATE(POINT_ZONE_INFO(MAX_ZONES,3))
      ALLOCATE(FACE_ZONE_INFO(MAX_ZONES,3))
      ALLOCATE(BC_LABEL_TEXT(MAX_ZONES,2))
      ALLOCATE(BC_ASSIGNED(MAX_ZONES))

      WRITE(ERR_MSG,2000) 'READING geometry from geometry.msh...'
      CALL LOG_INFO()

      INQUIRE(FILE='geometry.msh',EXIST=PRESENT)
      IF(.NOT.PRESENT) THEN
         IF(MyPE == PE_IO) THEN
            WRITE(ERR_MSG,"('(PE ',I3,'): input data file, ', A, ' is missing: run aborted')") &
            myPE,'geometry.msh'
         ENDIF
         CALL LOG_ERROR()
      ENDIF

! OPEN geometry.msh ASCII FILE
      OPEN(UNIT=333, FILE='geometry.msh', STATUS='OLD', ERR=910)
      WRITE(ERR_MSG,2000)'MSH file opened. Starting reading data...'
      CALL LOG_INFO()

      OPEN(UNIT=444, FILE='checkgeometry.stl')
      write(444,*)'solid vcg'

      CALL SKIP(333,3)

! Reading grid dimension
      READ(333,*) (BUFF_CHAR(I),I=1,2)
      CALL TEXT_HEX2INT(BUFF_CHAR(2)(1:1),GRID_DIMENSION)
      CALL SKIP(333,1)

! Reading Points
      READ(333,*) (BUFF_CHAR(I),I=1,4)
      CALL TEXT_HEX2INT(BUFF_CHAR(4),N_POINTS)
      ZONE = 1
      ALL_POINTS_READ = .FALSE.

      DO WHILE(.NOT.ALL_POINTS_READ)
         READ(333,*) (BUFF_CHAR(I),I=1,4)
         CALL TEXT_HEX2INT(BUFF_CHAR(2),ZONE_ID)
         CALL TEXT_HEX2INT(BUFF_CHAR(3),N1)
         CALL TEXT_HEX2INT(BUFF_CHAR(4),N2)
         POINT_ZONE_INFO(ZONE,1) = ZONE_ID
         POINT_ZONE_INFO(ZONE,2) = N1
         POINT_ZONE_INFO(ZONE,3) = N2

         DO NN = N1,N2
            READ(333,*) (POINT_COORDS(D,NN),D=1,GRID_DIMENSION)
         ENDDO

         IF(N2/=N_POINTS) THEN
            ZONE = ZONE + 1
         ELSE
            ALL_POINTS_READ = .TRUE.
!            WRITE(*,*)' ALL POINTS READ SUCCESSFULLY.'
          ENDIF
      ENDDO

      CALL SKIP(333,3)

! Reading FACES
      READ(333,*) (BUFF_CHAR(I),I=1,4)
      CALL TEXT_HEX2INT(BUFF_CHAR(3),N_FACES)
!      WRITE(*,*)' READING FACES ...'
      ZONE = 1
      ALL_FACES_READ = .FALSE.
      BC_LABELS_TO_READ = 0
      N_BC_IGNORED = 0
      NF=0
      N_QUAD2TRI=0

      DO WHILE(.NOT.ALL_FACES_READ)

         READ(333,*) (BUFF_CHAR(I),I=1,4)
         CALL TEXT_HEX2INT(BUFF_CHAR(1)(5:),ZONE_ID)
         CALL TEXT_HEX2INT(BUFF_CHAR(2),N1)
         CALL TEXT_HEX2INT(BUFF_CHAR(3),N2)

         FACE_ZONE_INFO(ZONE,1) = ZONE_ID
         FACE_ZONE_INFO(ZONE,2) = N1
         FACE_ZONE_INFO(ZONE,3) = N2
         BC_LABELS_TO_READ = BC_LABELS_TO_READ + 1

         IF(IS_CG(BC_TYPE_ENUM(ZONE_ID))) THEN

            DO NN = N1,N2
               READ(333,*)PPFACE
               IF(PPFACE<3.OR.PPFACE>4) THEN
                  WRITE(ERR_MSG,*)PPFACE, 'POINTS PER FACE. EACH FACE MUST HAVE 3 OR 4 POINTS.'
                  CALL LOG_ERROR()

               ELSEIF(PPFACE==3) THEN
                  BACKSPACE(333)
                  READ(333,*) (BUFF_CHAR(I),I=1,4)
                  CALL TEXT_HEX2INT(BUFF_CHAR(2),P1)
                  CALL TEXT_HEX2INT(BUFF_CHAR(3),P2)
                  CALL TEXT_HEX2INT(BUFF_CHAR(4),P3)

                  V1x = POINT_COORDS(1,P1)
                  V1y = POINT_COORDS(2,P1)
                  V1z = POINT_COORDS(3,P1)
                  V2x = POINT_COORDS(1,P2)
                  V2y = POINT_COORDS(2,P2)
                  V2z = POINT_COORDS(3,P2)
                  V3x = POINT_COORDS(1,P3)
                  V3y = POINT_COORDS(2,P3)
                  V3z = POINT_COORDS(3,P3)

                  VECTMP  = POINT_COORDS(:,P2)-POINT_COORDS(:,P1)
                  VECTMP2 = POINT_COORDS(:,P3)-POINT_COORDS(:,P1)
                  NORMAL = CROSS_PRODUCT(VECTMP,VECTMP2)
                  NORM = sqrt(dot_product(normal(:),normal(:)))
                  NORMAL = NORMAL / NORM

                  NF = NF + 1

! Save and Reverse unit vector if needed (this will switch fluid and blocked cells)
                  NORM_FACE(:,NF) = NORMAL*OUT_MSH_VALUE

                  VERTEX(1,1,NF) = SCALE_MSH*V1x + TX_MSH
                  VERTEX(1,2,NF) = SCALE_MSH*V1y + TY_MSH
                  VERTEX(1,3,NF) = SCALE_MSH*V1z + TZ_MSH
                  VERTEX(2,1,NF) = SCALE_MSH*V2x + TX_MSH
                  VERTEX(2,2,NF) = SCALE_MSH*V2y + TY_MSH
                  VERTEX(2,3,NF) = SCALE_MSH*V2z + TZ_MSH
                  VERTEX(3,1,NF) = SCALE_MSH*V3x + TX_MSH
                  VERTEX(3,2,NF) = SCALE_MSH*V3y + TY_MSH
                  VERTEX(3,3,NF) = SCALE_MSH*V3z + TZ_MSH

                  BC_ID_STL_FACE(NF) = ZONE_ID

                  write(444,*) '   facet normal ', NORM_FACE(:,NF)
                  write(444,*) '      outer loop'
                  write(444,*) '         vertex ', VERTEX(1,1:3,NF)
                  write(444,*) '         vertex ', VERTEX(2,1:3,NF)
                  write(444,*) '         vertex ', VERTEX(3,1:3,NF)
                  write(444,*) '      endloop'
                  write(444,*) '   endfacet'

               ELSEIF(PPFACE==4) THEN
                  BACKSPACE(333)
                  READ(333,*) (BUFF_CHAR(I),I=1,5)
                  CALL TEXT_HEX2INT(BUFF_CHAR(2),P1)
                  CALL TEXT_HEX2INT(BUFF_CHAR(3),P2)
                  CALL TEXT_HEX2INT(BUFF_CHAR(4),P3)
                  CALL TEXT_HEX2INT(BUFF_CHAR(5),P4)

                  N_QUAD2TRI = N_QUAD2TRI + 1

! Splitting Quad face 1-2-3-4 into two triangles 1-2-3 and 1-3-4

! First triangle 1-2-3
                  V1x = POINT_COORDS(1,P1)
                  V1y = POINT_COORDS(2,P1)
                  V1z = POINT_COORDS(3,P1)
                  V2x = POINT_COORDS(1,P2)
                  V2y = POINT_COORDS(2,P2)
                  V2z = POINT_COORDS(3,P2)
                  V3x = POINT_COORDS(1,P3)
                  V3y = POINT_COORDS(2,P3)
                  V3z = POINT_COORDS(3,P3)

                  VECTMP  = POINT_COORDS(:,P2)-POINT_COORDS(:,P1)
                  VECTMP2 = POINT_COORDS(:,P3)-POINT_COORDS(:,P1)
                  NORMAL = CROSS_PRODUCT(VECTMP,VECTMP2)
                  NORM = sqrt(dot_product(normal(:),normal(:)))
                  NORMAL = NORMAL / NORM

                  NF = NF + 1

                  ! Save and Reverse unit vector if needed (this will switch fluid and blocked cells)
                  NORM_FACE(:,NF) = NORMAL*OUT_MSH_VALUE
                  VERTEX(1,1,NF) = SCALE_MSH*V1x + TX_MSH
                  VERTEX(1,2,NF) = SCALE_MSH*V1y + TY_MSH
                  VERTEX(1,3,NF) = SCALE_MSH*V1z + TZ_MSH
                  VERTEX(2,1,NF) = SCALE_MSH*V2x + TX_MSH
                  VERTEX(2,2,NF) = SCALE_MSH*V2y + TY_MSH
                  VERTEX(2,3,NF) = SCALE_MSH*V2z + TZ_MSH
                  VERTEX(3,1,NF) = SCALE_MSH*V3x + TX_MSH
                  VERTEX(3,2,NF) = SCALE_MSH*V3y + TY_MSH
                  VERTEX(3,3,NF) = SCALE_MSH*V3z + TZ_MSH

                  BC_ID_STL_FACE(NF) = ZONE_ID


                  write(444,*) '   facet normal ', NORM_FACE(:,NF)
                  write(444,*) '      outer loop'
                  write(444,*) '         vertex ', VERTEX(1,1:3,NF)
                  write(444,*) '         vertex ', VERTEX(2,1:3,NF)
                  write(444,*) '         vertex ', VERTEX(3,1:3,NF)
                  write(444,*) '      endloop'
                  write(444,*)'   endfacet'

! Second triangle 1-3-4
                  V1x = POINT_COORDS(1,P1)
                  V1y = POINT_COORDS(2,P1)
                  V1z = POINT_COORDS(3,P1)
                  V2x = POINT_COORDS(1,P3)
                  V2y = POINT_COORDS(2,P3)
                  V2z = POINT_COORDS(3,P3)
                  V3x = POINT_COORDS(1,P4)
                  V3y = POINT_COORDS(2,P4)
                  V3z = POINT_COORDS(3,P4)

                  VECTMP  = POINT_COORDS(:,P3)-POINT_COORDS(:,P1)
                  VECTMP2 = POINT_COORDS(:,P4)-POINT_COORDS(:,P1)
                  NORMAL = CROSS_PRODUCT(VECTMP,VECTMP2)
                  NORM = sqrt(dot_product(normal(:),normal(:)))
                  NORMAL = NORMAL / NORM

                  NF = NF + 1

                  ! Save and Reverse unit vector if needed (this will switch fluid and blocked cells)
                  NORM_FACE(:,NF) = NORMAL*OUT_MSH_VALUE

                  VERTEX(1,1,NF) = SCALE_MSH*V1x + TX_MSH
                  VERTEX(1,2,NF) = SCALE_MSH*V1y + TY_MSH
                  VERTEX(1,3,NF) = SCALE_MSH*V1z + TZ_MSH
                  VERTEX(2,1,NF) = SCALE_MSH*V2x + TX_MSH
                  VERTEX(2,2,NF) = SCALE_MSH*V2y + TY_MSH
                  VERTEX(2,3,NF) = SCALE_MSH*V2z + TZ_MSH
                  VERTEX(3,1,NF) = SCALE_MSH*V3x + TX_MSH
                  VERTEX(3,2,NF) = SCALE_MSH*V3y + TY_MSH
                  VERTEX(3,3,NF) = SCALE_MSH*V3z + TZ_MSH

                  BC_ID_STL_FACE(NF) = ZONE_ID

                  write(444,*) '   facet normal ', NORM_FACE(:,NF)
                  write(444,*) '      outer loop'
                  write(444,*) '         vertex ', VERTEX(1,1:3,NF)
                  write(444,*) '         vertex ', VERTEX(2,1:3,NF)
                  write(444,*) '         vertex ', VERTEX(3,1:3,NF)
                  write(444,*) '      endloop'
                  write(444,*)'   endfacet'

               ENDIF
            ENDDO
         ELSE
            N_BC_IGNORED = N_BC_IGNORED + 1
!            WRITE(*,*)'BOUNDARY CONDITION FOR ZONE',ZONE_ID,' IS NOT DEFINED.'
!            WRITE(*,*)'THIS ZONE IS IGNORED'
            CALL SKIP(333,N2-N1+1)
         ENDIF


         IF(N2/=N_FACES) THEN
            ZONE = ZONE + 1
            ALL_FACES_READ = .FALSE.
            CALL SKIP(333,1)
         ELSE
            ALL_FACES_READ = .TRUE.
!            WRITE(*,*)' ALL FACES READ SUCCESSFULLY.'
            N_FACE_ZONES = ZONE
          ENDIF
      ENDDO

      write(444,*)'endsolid vcg'
      close(444)

      WRITE(ERR_MSG,*) '' // nl // &
         ' The file check_geometry.stl was successfully written.' // nl // &
         ' This is the equivalent of geometry.msh in STL format,' // nl // &
         ' and is provided for convenience (it is not used).'
      CALL LOG_INFO()

! Reading Boundary condition labels
      BC_LABELS_READ = 0
      BC_ASSIGNED = .FALSE.

      DO WHILE(BC_LABELS_READ < BC_LABELS_TO_READ)
         READ(333,*) BUFF_CHAR(1)
         IF(BUFF_CHAR(1)=='(45') THEN
            BACKSPACE(333)
            READ(333,*) (BUFF_CHAR(I),I=1,4)
            CALL TEXT_DEC2INT(BUFF_CHAR(2),ZONE_ID)
            DO ZONE = 1,N_FACE_ZONES
               IF(FACE_ZONE_INFO(ZONE,1)==ZONE_ID) THEN
                  BC_LABELS_READ = BC_LABELS_READ + 1
                  BC_LABEL_TEXT(ZONE,1) = TRIM(BUFF_CHAR(3))
                  BC_LABEL_TEXT(ZONE,2) = TRIM(BUFF_CHAR(4))
                  IF(IS_CG(BC_TYPE_ENUM(ZONE_ID))) THEN
                     BC_ASSIGNED(ZONE) = .TRUE.
                  ENDIF
               ENDIF
            ENDDO
         ENDIF
      ENDDO

      IF(MyPE == PE_IO) THEN
         WRITE(*,*)' Summary of data read from geometry.msh file:'
         WRITE(*,*)'======================================================================'
         WRITE(*,*)' DIMENSION   POINTS       FACES       ZONES'
         WRITE(*,*)'======================================================================'
         WRITE(*,1020) GRID_DIMENSION,N_POINTS,N_FACES,N_FACE_ZONES
         WRITE(*,*)''
         WRITE(*,*)' BOUNDARY CONDITION DETECTED (INFO EXTRACTED FROM .MSH FILE):'
         WRITE(*,*)'======================================================================'
         WRITE(*,*)'  ZONE  BC_TYPE(MFIX)  FACES   BC_TYPE(GAMBIT)      BC LABEL'
         WRITE(*,*)'======================================================================'

         N_FACES = 0
         DO ZONE = 1,N_FACE_ZONES
            IF(BC_ASSIGNED(ZONE)) THEN
               ZID = FACE_ZONE_INFO(ZONE,1)
               NFZ = FACE_ZONE_INFO(ZONE,3)-FACE_ZONE_INFO(ZONE,2) + 1
               N_FACES = N_FACES + NFZ
               L2=LEN(TRIM(BC_LABEL_TEXT(ZONE,2)))-4
               WRITE(*,1000) ZID,BC_TYPE_ENUM(ZID),NFZ,BC_LABEL_TEXT(ZONE,1),BC_LABEL_TEXT(ZONE,2)(1:L2)
            ENDIF
         ENDDO

         WRITE(*,*)''
         DO ZONE = 1,N_FACE_ZONES
            IF(.NOT.BC_ASSIGNED(ZONE)) THEN
               ZID = FACE_ZONE_INFO(ZONE,1)
               NFZ = FACE_ZONE_INFO(ZONE,3)-FACE_ZONE_INFO(ZONE,2) + 1
               L2=LEN(TRIM(BC_LABEL_TEXT(ZONE,2)))-4
               WRITE(*,1000) ZID,'NOT USED',NFZ,BC_LABEL_TEXT(ZONE,1),BC_LABEL_TEXT(ZONE,2)(1:L2)
            ENDIF
         ENDDO
         WRITE(*,*)'======================================================================'
         WRITE(*,*)' PLEASE VERIFY THAT BOUNDARY CONDITIONS ARE CORRECTLY ASSIGNED.'
         WRITE(*,*)' MODIFY BC_TYPE keyword IF NECESSARY.'
         WRITE(*,*)''
      ENDIF   ! endif mype=pe_io

      XMIN_MSH = MINVAL(VERTEX(:,1,1:N_FACES))
      XMAX_MSH = MAXVAL(VERTEX(:,1,1:N_FACES))
      YMIN_MSH = MINVAL(VERTEX(:,2,1:N_FACES))
      YMAX_MSH = MAXVAL(VERTEX(:,2,1:N_FACES))
      ZMIN_MSH = MINVAL(VERTEX(:,3,1:N_FACES))
      ZMAX_MSH = MAXVAL(VERTEX(:,3,1:N_FACES))

      IF(MyPE == PE_IO) THEN
         WRITE(*,2000)'MSH file successfully read.'
         WRITE(*,*)' Total number of faces used as boundary faces =',N_FACES
         IF(N_QUAD2TRI>0) WRITE(*,*)' Number of quad faces split into triangles    =',N_QUAD2TRI


         WRITE(*,*)' RANGE OF MSH FILE:'
         IF(SCALE_MSH/=ONE) THEN
            WRITE(*,5000)' AFTER SCALING BY A FACTOR OF ',SCALE_MSH
         ENDIF
         ABSTRANS = dabs(TX_MSH)+dabs(TY_MSH)+dabs(TZ_MSH)
         IF(ABSTRANS>TOL_MSH) THEN
            WRITE(*,3000)' AFTER TRANSLATION OF (X,Y,Z)=',TX_MSH,TY_MSH,TZ_MSH
         ENDIF
         WRITE(*,4000)'X-RANGE = ', XMIN_MSH,XMAX_MSH
         WRITE(*,4000)'Y-RANGE = ', YMIN_MSH,YMAX_MSH
         WRITE(*,4000)'Z-RANGE = ', ZMIN_MSH,ZMAX_MSH
         WRITE(*,4000)''
      ENDIF

      XMIN_MSH = XMIN_MSH - 10.0*TOL_MSH
      XMAX_MSH = XMAX_MSH + 10.0*TOL_MSH
      YMIN_MSH = YMIN_MSH - 10.0*TOL_MSH
      YMAX_MSH = YMAX_MSH + 10.0*TOL_MSH
      ZMIN_MSH = ZMIN_MSH - 10.0*TOL_MSH
      ZMAX_MSH = ZMAX_MSH + 10.0*TOL_MSH

      N_FACETS = N_FACES

      N_FACETS = NF

      XMIN_STL = XMIN_MSH
      XMAX_STL = XMAX_MSH
      YMIN_STL = YMIN_MSH
      YMAX_STL = YMAX_MSH
      ZMIN_STL = ZMIN_MSH
      ZMAX_STL = ZMAX_MSH

      OUT_STL_VALUE = OUT_MSH_VALUE

      CLOSE(333)

      OPEN(UNIT=444, FILE='msgeometry.stl')

         DO ZONE = 1,N_FACE_ZONES
            IF(BC_ASSIGNED(ZONE)) THEN
               ZID = FACE_ZONE_INFO(ZONE,1)
               L2=LEN(TRIM(BC_LABEL_TEXT(ZONE,2)))-4
!               WRITE(*,1000) ZID,BC_TYPE(ZID),NFZ,BC_LABEL_TEXT(ZONE,1),BC_LABEL_TEXT(ZONE,2)(1:L2)
               print*,'=======>  Zone ID= ', ZID,BC_LABEL_TEXT(ZONE,2)(1:L2)
               write(444,6000) 'solid', BC_LABEL_TEXT(ZONE,2)(1:L2),ZID
               DO NF=1,N_FACETS
                  IF(BC_ID_STL_FACE(NF) == ZID) THEN

                  write(444,*) '   facet normal ', NORM_FACE(:,NF)
                  write(444,*) '      outer loop'
                  write(444,*) '         vertex ', VERTEX(1,1:3,NF)
                  write(444,*) '         vertex ', VERTEX(2,1:3,NF)
                  write(444,*) '         vertex ', VERTEX(3,1:3,NF)
                  write(444,*) '      endloop'
                  write(444,*)'   endfacet'


                  ENDIF
               ENDDO

               write(444,6000) 'endsolid', BC_LABEL_TEXT(ZONE,2)(1:L2),ZID

            ENDIF
         ENDDO

      close(444)

      DEALLOCATE(POINT_COORDS)
      DEALLOCATE(POINT_ZONE_INFO)
      DEALLOCATE(FACE_ZONE_INFO)
      DEALLOCATE(BC_LABEL_TEXT)
      DEALLOCATE(BC_ASSIGNED)
      RETURN


! HERE IF AN ERROR OCCURRED OPENING/READING THE FILE
 910  CONTINUE
      WRITE (*, 1500)
      call log_error()
 920  CONTINUE
      WRITE (*, 1600)
      call log_error()
 930  CONTINUE
      WRITE (*, 1700)
      call log_error()

 1000 FORMAT(1X,I4,7X,A8,I8,1X,A20,1X,A20)
 1020 FORMAT(5X,I3,2X,3(I10,2X))
!
 1500 FORMAT(/1X,70('*')//' From: GET_MSH_DATA',/' Message: ',&
      'Unable to open msh file',/1X,70('*')/)
 1600 FORMAT(/1X,70('*')//' From: GET_MSH_DATA',/' Message: ',&
      'Error while reading msh file',/1X,70('*')/)
 1700 FORMAT(/1X,70('*')//' From: GET_MSH_DATA',/' Message: ',&
      'End of file reached while reading msh file',/1X,70('*')/)
 2000 FORMAT(1X,A)
 3000 FORMAT(1X,A,'(',F10.4,';',F10.4,';',F10.4,')')
 4000 FORMAT(1X,A,F10.4,' to ',F10.4)
 5000 FORMAT(1X,A,F10.4)
 6000 FORMAT(A,1X,A,'_',I4.4)
      END SUBROUTINE GET_MSH_DATA


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE SKIP(FILE_UNIT,N_SKIP)

      IMPLICIT NONE
      INTEGER, INTENT(IN) ::FILE_UNIT,N_SKIP
      INTEGER :: I
!---------------------------------------------------------------------

      DO I = 1,N_SKIP
         READ(FILE_UNIT,*)
      ENDDO
      RETURN
      END SUBROUTINE SKIP


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE TEXT_HEX2INT(STRING,INT)

      CHARACTER(LEN=10)  :: INTERNAL
      CHARACTER(LEN=*)   :: STRING
      CHARACTER(LEN=10)  :: CLEANED_STRING
      INTEGER :: INT
!---------------------------------------------------------------------

      CLEANED_STRING = STRING
      IF(STRING(1:1)=='(') CLEANED_STRING = STRING(2:)
      WRITE(INTERNAL,fmt='(A10)') TRIM(CLEANED_STRING)
      READ(INTERNAL,FMT='(Z10)') INT
      RETURN
      END SUBROUTINE TEXT_HEX2INT


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE TEXT_DEC2INT(STRING,INT)

      CHARACTER(LEN=10)  :: INTERNAL
      CHARACTER(LEN=*)   :: STRING
      CHARACTER(LEN=10)  :: CLEANED_STRING
      INTEGER :: INT
!---------------------------------------------------------------------

      CLEANED_STRING = STRING
      IF(STRING(1:1)=='(') CLEANED_STRING = STRING(2:)
      WRITE(INTERNAL,fmt='(A10)') TRIM(CLEANED_STRING)
      READ(INTERNAL,FMT='(I10)') INT
      RETURN
      END SUBROUTINE TEXT_DEC2INT


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: GET_STL_DATA                                           C
!  Purpose: reads face vertices  and normal vectors from an STL file   C
!                                                                      C
!  Author: Jeff Dietiker                              Date: 30-JAN-09  C
!  Reviewer:                                          Date: **-***-**  C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE GET_STL_DATA

      IMPLICIT NONE
      INTEGER :: NN, NG, IGNORED_FACETS
      LOGICAL :: PRESENT,KEEP_READING,IGNORE_CURRENT_FACET
      DOUBLE PRECISION ::v1x,v1y,v1z
      DOUBLE PRECISION ::v2x,v2y,v2z
      DOUBLE PRECISION ::v3x,v3y,v3z
      DOUBLE PRECISION ::x12,y12,z12,x13,y13,z13,dp,d12,d13
      DOUBLE PRECISION ::cos_angle,cos_small_angle
      DOUBLE PRECISION ::n1,n2,n3,NORM
      DOUBLE PRECISION ::ABSTRANS
      CHARACTER(LEN=32) ::TEST_CHAR,BUFF_CHAR
      CHARACTER(LEN=255) ::geometryfile(0:DIMENSION_BC)
      INTEGER :: BCV,NUMBER_OF_GEOMETRY_FILES,NUMBER_OF_BC_PATCHES
      INTEGER :: BC_PATCH(0:DIMENSION_BC),L2,NF1,NF2
      LOGICAL :: BC_PATCH_FOUND_IN_STL(DIMENSION_BC)
      CHARACTER(LEN=100) :: FNAME
      integer :: stl_unit, nf
      LOGICAL :: MULTISTL
!---------------------------------------------------------------------

      IF(.NOT.USE_STL) RETURN

      IF(NO_K) THEN
         WRITE(ERR_MSG,*) &
            'ERROR: STL METHOD VALID ONLY IN 3D.'
         CALL LOG_ERROR()
      ENDIF
      IF(N_QUADRIC > 0) THEN
         WRITE(ERR_MSG,*) 'ERROR: BOTH QUADRIC(S) AND '//&
            'STL INPUT ARE SPECIFIED.'
         WRITE(ERR_MSG,*) 'MFIX HANDLES ONLY ONE TYPE '//&
            'OF SURFACE INPUT.'
         CALL LOG_ERROR()
      ENDIF
      IF(USE_MSH) THEN
         WRITE(ERR_MSG,*) 'ERROR: BOTH MSH AND '//&
            'STL INPUT ARE SPECIFIED.'
         WRITE(ERR_MSG,*) 'MFIX HANDLES ONLY ONE TYPE '//&
            'OF SURFACE INPUT.'
         CALL LOG_ERROR()
      ENDIF

      MULTISTL = .FALSE.
      BC_PATCH_FOUND_IN_STL = .FALSE.
      NUMBER_OF_GEOMETRY_FILES = 0
      NUMBER_OF_BC_PATCHES     = 0
      N_STL_GROUP = 0

! DETERMINE WHICH BOUNDARY CONDITIONS NEED STL FILE
      WRITE(ERR_MSG,100)'From Get_Stl_Data: Analyzing BCs...'
      CALL LOG_STATUS()
      WRITE(ERR_MSG,120)'BC_ID','BC_TYPE'
      CALL LOG_STATUS()

      DO BCV = 1, DIMENSION_BC
!         IF(IS_CG(BC_TYPE_ENUM(BCV))) THEN
         IF(BC_TYPE(BCV)(1:3)=='CG_') THEN
            NUMBER_OF_GEOMETRY_FILES = NUMBER_OF_GEOMETRY_FILES + 1
            NUMBER_OF_BC_PATCHES     = NUMBER_OF_BC_PATCHES + 1
            BC_PATCH(NUMBER_OF_GEOMETRY_FILES) = BCV
            WRITE(geometryfile(NUMBER_OF_GEOMETRY_FILES),200) 'geometry_',BCV
            INQUIRE(FILE=TRIM(geometryfile(NUMBER_OF_GEOMETRY_FILES)),&
               EXIST=PRESENT)
            MULTISTL = MULTISTL .OR. PRESENT
            WRITE(ERR_MSG,130)BCV,BC_TYPE(BCV)
            CALL LOG_STATUS()
         ENDIF
      ENDDO

      WRITE(ERR_MSG,110)'Number of CG BCs = ',NUMBER_OF_BC_PATCHES
      CALL LOG_STATUS()

100  FORMAT(1X,A)
110  FORMAT(1X,A,I6)
120  FORMAT(1X,A6,4X,A7)
130  FORMAT(1X,I6,4X,A6)
140  FORMAT(1X,A,I6,A)
200  FORMAT(A,I4.4,".stl")

!!!  TODO cleanup this message cgw

! VERIFY THAT THERE IS AT LEAST ONE STL FILE TO READ
      IF(NUMBER_OF_BC_PATCHES==0) THEN
         IF(MyPE == PE_IO) THEN
            WRITE(ERR_MSG,'(A)') 'ERROR: NO CARTESIAN GRID BOUNDARY CONDITION SPECIFIED.' // new_line('') // &
               'AT LEAST ONE BC_TYPE MUST START WITH CG (FOR EXAMPLE CG_NSW)' // new_line('') // &
               'RUN ABORTED'
            CALL LOG_ERROR()
         ENDIF
      ENDIF


! START READING EACH STL FILE, ONE FOR EACH BC_TYPE
      N_FACETS = 0
      IGNORED_FACETS = 0

! Open file that stores bad facets (they are ignored)
      IF(mype.eq.pe_io) then
         WRITE(fname,'(A,"_FACETS_IGNORED", ".stl")') &
         TRIM(RUN_NAME)
         open(newunit=stl_unit, file = fname, form='formatted')
         write(stl_unit,*)'solid vcg'
      ENDIF


! Read each geometry_####.stl file
      DO NN = 1, NUMBER_OF_GEOMETRY_FILES
         WRITE(ERR_MSG,2000) 'Reading geometry from '//TRIM(geometryfile(NN))//' ...'
         CALL LOG_INFO()
         INQUIRE(FILE=TRIM(geometryfile(NN)),EXIST=PRESENT)
         IF(.NOT.PRESENT) THEN
            WRITE(ERR_MSG,"('(PE ',I3,'): input data file, ',A20,' is missing: run aborted')") &
               myPE,TRIM(geometryfile(NN))
            CALL LOG_ERROR()
         ENDIF

         N_STL_GROUP                     = N_STL_GROUP + 1
         STL_GROUP(N_STL_GROUP)%TYPE     = 'BC_STL'
         STL_GROUP(N_STL_GROUP)%FILENAME = geometryfile(NN)
         STL_GROUP(N_STL_GROUP)%ID       = BC_PATCH(NN)
         IF(N_STL_GROUP==1) THEN
            STL_GROUP(N_STL_GROUP)%START    = 1
         ELSE
            STL_GROUP(N_STL_GROUP)%START    = STL_GROUP(N_STL_GROUP - 1)%END + 1
         ENDIF

! OPEN geometry.stl ASCII FILE
         OPEN(UNIT=333, FILE=TRIM(geometryfile(NN)), STATUS='OLD', ERR=910)
         WRITE(ERR_MSG,2000)'STL file opened. Starting reading data...'
         CALL LOG_STATUS()
         KEEP_READING = .TRUE.
         DO WHILE(KEEP_READING)

            READ(333,*,ERR=920,END=930) TEST_CHAR
!            print *,'TEST_CHAR=',TEST_CHAR
            IF(TRIM(TEST_CHAR) == 'solid'.AND.NN==0) THEN

               BACKSPACE(333)
               READ(333,*,ERR=920,END=930) BUFF_CHAR,BUFF_CHAR
               L2=LEN(TRIM(BUFF_CHAR))-3
               READ(BUFF_CHAR(L2:L2+4),fmt='(I4.4)') BC_PATCH(NN)
               WRITE(ERR_MSG,140) 'Found BC patch #', BC_PATCH(NN), ' in STL file.'
               call log_info()
               IF(.NOT.IS_CG(BC_TYPE_ENUM(BC_PATCH(NN)))) THEN
                  IF(MyPE == PE_IO)  THEN
                     WRITE(ERR_MSG,110) 'This BC patch does not match a CG BC:',BC_PATCH(NN)
                     call log_error()
                     return
                  ENDIF
               ENDIF
               BC_PATCH_FOUND_IN_STL(BC_PATCH(NN)) = .TRUE.

            ELSEIF(TRIM(TEST_CHAR) == 'facet') THEN

               BACKSPACE(333)
               IGNORE_CURRENT_FACET = .FALSE.
               READ(333,*,ERR=920,END=930) BUFF_CHAR,BUFF_CHAR,N1,N2,N3  ! Read unit normal vector
               READ(333,*,ERR=920,END=930)
               READ(333,*,ERR=920,END=930) BUFF_CHAR, V1x,V1y,V1z
               READ(333,*,ERR=920,END=930) BUFF_CHAR, V2x,V2y,V2z
               READ(333,*,ERR=920,END=930) BUFF_CHAR, V3x,V3y,V3z
               N1 = N1 * OUT_STL_VALUE  ! Reverse unit vector if needed (this will switch fluid and blocked cells)
               N2 = N2 * OUT_STL_VALUE
               N3 = N3 * OUT_STL_VALUE

               IF(FLIP_STL_NORMALS(BC_PATCH(NN))) THEN
                  N1 = - N1
                  N2 = - N2
                  N3 = - N3
               ENDIF

               x12 = V2x - V1x
               y12 = V2y - V1y
               z12 = V2z - V1z
               x13 = V3x - V1x
               y13 = V3y - V1y
               z13 = V3z - V1z
               dp  = x12*x13 + y12*y13 + z12*z13
               d12 = sqrt(x12**2+y12**2+z12**2)
               d13 = sqrt(x13**2+y13**2+z13**2)

               IF((d12*d13)>TOL_STL**2) THEN
                  cos_angle = dp/(d12*d13)
               ELSE
                  cos_angle = ONE
               ENDIF
               cos_small_angle = dcos(STL_SMALL_ANGLE / 180.0 * PI)

               IF(DABS(cos_angle)>cos_small_angle) THEN
                  IGNORE_CURRENT_FACET = .TRUE.    ! Ignore small facets
               ENDIF

               NORM = sqrt(N1**2+N2**2+N3**2)

               IF(NORM>TOL_STL) THEN
                  N1 = N1 /NORM
                  N2 = N2 /NORM
                  N3 = N3 /NORM
               ELSE
                  IGNORE_CURRENT_FACET = .TRUE.  ! Ignore facets with zero normal vector
               ENDIF

               IF(IGNORE_CURRENT_FACET) THEN
                  IGNORED_FACETS = IGNORED_FACETS + 1
                  IF(mype.eq.pe_io) then
                     write(stl_unit,*) '   facet normal ', N1, N2, N3
                     write(stl_unit,*) '      outer loop'
                     write(stl_unit,*) '         vertex ', V1x,V1y,V1z
                     write(stl_unit,*) '         vertex ', V2x,V2y,V2z
                     write(stl_unit,*) '         vertex ', V3x,V3y,V3z
                     write(stl_unit,*) '      endloop'
                     write(stl_unit,*)'   endfacet'
                  ENDIF

               ELSE                              ! Save vertex coordinates for valid facets
                  N_FACETS = N_FACETS + 1
                  STL_GROUP(N_STL_GROUP)%END = N_FACETS
                  NORM_FACE(1,N_FACETS) = N1
                  NORM_FACE(2,N_FACETS) = N2
                  NORM_FACE(3,N_FACETS) = N3
                  VERTEX(1,1,N_FACETS) = SCALE_STL*V1x + TX_STL
                  VERTEX(1,2,N_FACETS) = SCALE_STL*V1y + TY_STL
                  VERTEX(1,3,N_FACETS) = SCALE_STL*V1z + TZ_STL
                  VERTEX(2,1,N_FACETS) = SCALE_STL*V2x + TX_STL
                  VERTEX(2,2,N_FACETS) = SCALE_STL*V2y + TY_STL
                  VERTEX(2,3,N_FACETS) = SCALE_STL*V2z + TZ_STL
                  VERTEX(3,1,N_FACETS) = SCALE_STL*V3x + TX_STL
                  VERTEX(3,2,N_FACETS) = SCALE_STL*V3y + TY_STL
                  VERTEX(3,3,N_FACETS) = SCALE_STL*V3z + TZ_STL
                  BC_ID_STL_FACE(N_FACETS) = BC_PATCH(NN)

               ENDIF

!            ELSEIF(TRIM(TEST_CHAR) == 'endsolid') THEN
!               KEEP_READING = .FALSE.
            ENDIF
         ENDDO

930      WRITE(ERR_MSG,100) 'Done reading file.'
         CALL LOG_STATUS()
         CLOSE(333)
      ENDDO


      XMIN_STL = MINVAL(VERTEX(:,1,1:N_FACETS))
      XMAX_STL = MAXVAL(VERTEX(:,1,1:N_FACETS))
      YMIN_STL = MINVAL(VERTEX(:,2,1:N_FACETS))
      YMAX_STL = MAXVAL(VERTEX(:,2,1:N_FACETS))
      ZMIN_STL = MINVAL(VERTEX(:,3,1:N_FACETS))
      ZMAX_STL = MAXVAL(VERTEX(:,3,1:N_FACETS))

      WRITE(ERR_MSG,2000)'STL file(s) successfully read.'
      CALL LOG_STATUS()
      WRITE(ERR_MSG,110)'Total number of facets read =',N_FACETS + IGNORED_FACETS
      CALL LOG_STATUS()
      WRITE(ERR_MSG,110)'Number of valid facets      =',N_FACETS
      CALL LOG_STATUS()
      WRITE(ERR_MSG,110)'Number of ignored facets    =',IGNORED_FACETS
      CALL LOG_STATUS()
      WRITE(ERR_MSG,100)'Geometry range from stl file:'
      CALL LOG_STATUS()
      IF(SCALE_STL/=ONE) THEN
         WRITE(ERR_MSG,5000)'Scaling factor:',SCALE_STL
      ENDIF
      ABSTRANS = dabs(TX_STL)+dabs(TY_STL)+dabs(TZ_STL)
      IF(ABSTRANS>TOL_STL) THEN
         WRITE(ERR_MSG,3000)'Translation vector (X,Y,Z)=',TX_STL,TY_STL,TZ_STL
         CALL LOG_STATUS()
      ENDIF
      WRITE(ERR_MSG,4000)'x-range = ', XMIN_STL,XMAX_STL
      CALL LOG_STATUS()
      WRITE(ERR_MSG,4000)'y-range = ', YMIN_STL,YMAX_STL
      CALL LOG_STATUS()
      WRITE(ERR_MSG,4000)'z-range = ', ZMIN_STL,ZMAX_STL
      CALL LOG_STATUS()
      WRITE(ERR_MSG,4000)''
      CALL LOG_STATUS()

      XMIN_STL = XMIN_STL - 10.0*TOL_STL
      XMAX_STL = XMAX_STL + 10.0*TOL_STL
      YMIN_STL = YMIN_STL - 10.0*TOL_STL
      YMAX_STL = YMAX_STL + 10.0*TOL_STL
      ZMIN_STL = ZMIN_STL - 10.0*TOL_STL
      ZMAX_STL = ZMAX_STL + 10.0*TOL_STL

!      IF(XMIN_STL<ZERO) XMIN_STL=ZERO
!      IF(XMAX_STL>XLENGTH) XMAX_STL=XLENGTH
!      IF(YMIN_STL<ZERO) YMIN_STL=ZERO
!      IF(YMAX_STL>YLENGTH) YMAX_STL=YLENGTH
!      IF(ZMIN_STL<ZERO) ZMIN_STL=ZERO
!      IF(ZMAX_STL>ZLENGTH) ZMAX_STL=ZLENGTH

! Close FACETS_IGNORED.stl
      IF(mype.eq.pe_io) then
         write(stl_unit,*)'endsolid vcg'
         close(stl_unit, status = 'keep')
         WRITE(ERR_MSG,100) 'The file FACETS_IGNORED.stl was successfully written.'
         CALL LOG_STATUS()
         WRITE(ERR_MSG,100) 'and is provided for convenience (it is not used).'
         CALL LOG_STATUS()
      ENDIF

      CALL DETECT_STL_ALONG_MFIX_BOX(NUMBER_OF_BC_PATCHES)

      IF(mype.eq.pe_io.and..true.) then
         WRITE(fname,'(A,"_FACETS_READ", ".stl")') &
         TRIM(RUN_NAME)
         open(newunit=stl_unit, file = fname, form='formatted')
         write(stl_unit,*)'solid vcg'

         do nf = 1, n_facets
            if(.not.stl_facet_along_box(nf)) then
               write(stl_unit,*) '   facet normal ', NORM_FACE(:,NF)
               write(stl_unit,*) '      outer loop'
               write(stl_unit,*) '         vertex ', VERTEX(1,1:3,NF)
               write(stl_unit,*) '         vertex ', VERTEX(2,1:3,NF)
               write(stl_unit,*) '         vertex ', VERTEX(3,1:3,NF)
               write(stl_unit,*) '      endloop'
               write(stl_unit,*)'   endfacet'
            endif
         enddo

         write(stl_unit,*)'endsolid vcg'
         close(stl_unit, status = 'keep')
         WRITE(ERR_MSG,100) 'The file FACETS_READ.stl was successfully written.'
         CALL LOG_STATUS()
         WRITE(ERR_MSG,100) 'and is provided for convenience (it is not used).'
         CALL LOG_STATUS()
      ENDIF


      IF(.FALSE.) THEN
         DO NG=1,N_STL_GROUP
            WRITE(*,*) 'STL GROUP: ', NG
            WRITE(*,*) 'TYPE     : ', STL_GROUP(NG)%TYPE
            WRITE(*,*) 'FILENAME : ', STL_GROUP(NG)%FILENAME
            WRITE(*,*) 'ID       : ', STL_GROUP(NG)%ID
            WRITE(*,*) 'START    : ', STL_GROUP(NG)%START
            WRITE(*,*) 'END      : ', STL_GROUP(NG)%END
         ENDDO
      ENDIF


      RETURN
! HERE IF AN ERROR OCCURRED OPENING/READING THE FILE

 910  CONTINUE
      WRITE (ERR_MSG, 1500)
      call log_error()
 920  CONTINUE
      WRITE (ERR_MSG, 1600)
      call log_error()
! 930  CONTINUE
!     WRITE (ERR_MSG, 1700)
!     call log_error()

 1500 FORMAT(/1X,70('*')//' From: GET_STL_DATA',/' Message: ',&
      'Unable to open stl file',/1X,70('*')/)
 1600 FORMAT(/1X,70('*')//' From: GET_STL_DATA',/' Message: ',&
      'Error while reading stl file',/1X,70('*')/)
 2000 FORMAT(1X,A)
 3000 FORMAT(1X,A,'(',F10.4,';',F10.4,';',F10.4,')')
 4000 FORMAT(1X,A,F10.4,' to ',F10.4)
 5000 FORMAT(1X,A,F10.4)
      END SUBROUTINE GET_STL_DATA

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: GET_STL_DATA_IS                                        C
!  Purpose: reads face vertices  and normal vectors from an STL file   C
!           This version reads internal surfaces stl files only        C
!                                                                      C
!  Author: Jeff Dietiker                              Date: 09-JUN-20  C
!  Reviewer:                                          Date: **-***-**  C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE GET_STL_DATA_IS

      IMPLICIT NONE
      INTEGER :: NN, NG, IGNORED_FACETS
      LOGICAL :: PRESENT,KEEP_READING,IGNORE_CURRENT_FACET
      DOUBLE PRECISION ::v1x,v1y,v1z
      DOUBLE PRECISION ::v2x,v2y,v2z
      DOUBLE PRECISION ::v3x,v3y,v3z
      DOUBLE PRECISION ::x12,y12,z12,x13,y13,z13,dp,d12,d13
      DOUBLE PRECISION ::cos_angle,cos_small_angle
      DOUBLE PRECISION ::n1,n2,n3,NORM
      DOUBLE PRECISION ::ABSTRANS
      CHARACTER(LEN=32) ::TEST_CHAR,BUFF_CHAR
      CHARACTER(LEN=255) ::geometryfile(0:DIMENSION_IS)
      INTEGER :: ISV,NUMBER_OF_GEOMETRY_FILES,NUMBER_OF_IS_PATCHES
      INTEGER :: IS_PATCH(0:DIMENSION_IS),L2,NF1,NF2
      LOGICAL :: IS_PATCH_FOUND_IN_STL(DIMENSION_IS)
      CHARACTER(LEN=100) :: FNAME
      integer :: stl_unit, nf
      LOGICAL :: MULTISTL
!---------------------------------------------------------------------

      IF(.NOT.USE_STL) RETURN

      IF(NO_K) THEN
         WRITE(ERR_MSG,*) &
            'ERROR: STL METHOD VALID ONLY IN 3D.'
         CALL LOG_ERROR()
      ENDIF
      IF(N_QUADRIC > 0) THEN
         WRITE(ERR_MSG,*) 'ERROR: BOTH QUADRIC(S) AND '//&
            'STL INPUT ARE SPECIFIED.'
         WRITE(ERR_MSG,*) 'MFIX HANDLES ONLY ONE TYPE '//&
            'OF SURFACE INPUT.'
         CALL LOG_ERROR()
      ENDIF
      IF(USE_MSH) THEN
         WRITE(ERR_MSG,*) 'ERROR: BOTH MSH AND '//&
            'STL INPUT ARE SPECIFIED.'
         WRITE(ERR_MSG,*) 'MFIX HANDLES ONLY ONE TYPE '//&
            'OF SURFACE INPUT.'
         CALL LOG_ERROR()
      ENDIF

      NUMBER_OF_GEOMETRY_FILES = 0

! DETERMINE WHICH BOUNDARY CONDITIONS NEED STL FILE
      WRITE(ERR_MSG,100)'From Get_Stl_Data_is: Analyzing ISs...'
      CALL LOG_STATUS()
      WRITE(ERR_MSG,120)'IS_ID','IS_TYPE'
      CALL LOG_STATUS()

      DO ISV = 1, DIMENSION_IS
         IF(IS_TYPE(ISV)(1:3)=='STL') THEN
            NUMBER_OF_GEOMETRY_FILES = NUMBER_OF_GEOMETRY_FILES + 1
            IS_PATCH(NUMBER_OF_GEOMETRY_FILES) = ISV
            WRITE(geometryfile(NUMBER_OF_GEOMETRY_FILES),200) 'is_',ISV
            INQUIRE(FILE=TRIM(geometryfile(NUMBER_OF_GEOMETRY_FILES)),&
               EXIST=PRESENT)
            WRITE(ERR_MSG,130)ISV,IS_TYPE(ISV)
            CALL LOG_STATUS()
         ENDIF
      ENDDO

      WRITE(ERR_MSG,110)'Number of IS using STL = ',NUMBER_OF_GEOMETRY_FILES
      CALL LOG_STATUS()

100  FORMAT(1X,A)
110  FORMAT(1X,A,I6)
120  FORMAT(1X,A6,4X,A7)
130  FORMAT(1X,I6,4X,A6)
140  FORMAT(1X,A,I6,A)
200  FORMAT(A,I4.4,".stl")

      NF1 = 1
      NF2 = NUMBER_OF_GEOMETRY_FILES

! START READING EACH STL FILE, ONE FOR EACH IS_TYPE
      IGNORED_FACETS = 0

      DO NN = NF1, NF2
         WRITE(ERR_MSG,2000) 'Reading IS geometry from '//TRIM(geometryfile(NN))//' ...'
         CALL LOG_INFO()
         INQUIRE(FILE=TRIM(geometryfile(NN)),EXIST=PRESENT)
         IF(.NOT.PRESENT) THEN
            WRITE(ERR_MSG,"('(PE ',I3,'): input data file, ',A20,' is missing: run aborted')") &
               myPE,TRIM(geometryfile(NN))
            CALL LOG_ERROR()
         ENDIF

! Add data to STL_GROUP
         N_STL_GROUP                     = N_STL_GROUP + 1
         STL_GROUP(N_STL_GROUP)%TYPE     = 'IS_STL'
         STL_GROUP(N_STL_GROUP)%FILENAME = geometryfile(NN)
         STL_GROUP(N_STL_GROUP)%ID       = IS_PATCH(NN)
         IF(N_STL_GROUP==1) THEN
            STL_GROUP(N_STL_GROUP)%START    = 1
         ELSE
            STL_GROUP(N_STL_GROUP)%START    = STL_GROUP(N_STL_GROUP - 1)%END + 1
         ENDIF

! OPEN geometry.stl ASCII FILE
         OPEN(UNIT=333, FILE=TRIM(geometryfile(NN)), STATUS='OLD', ERR=910)
         WRITE(ERR_MSG,2000)'STL file opened. Starting reading data...'
         CALL LOG_STATUS()
         KEEP_READING = .TRUE.
         DO WHILE(KEEP_READING)

            READ(333,*,ERR=920,END=930) TEST_CHAR
            IF(TRIM(TEST_CHAR) == 'facet') THEN

               BACKSPACE(333)
               IGNORE_CURRENT_FACET = .FALSE.
               READ(333,*,ERR=920,END=930) BUFF_CHAR,BUFF_CHAR,N1,N2,N3  ! Read unit normal vector
               READ(333,*,ERR=920,END=930)
               READ(333,*,ERR=920,END=930) BUFF_CHAR, V1x,V1y,V1z
               READ(333,*,ERR=920,END=930) BUFF_CHAR, V2x,V2y,V2z
               READ(333,*,ERR=920,END=930) BUFF_CHAR, V3x,V3y,V3z

               x12 = V2x - V1x
               y12 = V2y - V1y
               z12 = V2z - V1z
               x13 = V3x - V1x
               y13 = V3y - V1y
               z13 = V3z - V1z
               dp  = x12*x13 + y12*y13 + z12*z13
               d12 = sqrt(x12**2+y12**2+z12**2)
               d13 = sqrt(x13**2+y13**2+z13**2)

               IF((d12*d13)>TOL_STL**2) THEN
                  cos_angle = dp/(d12*d13)
               ELSE
                  cos_angle = ONE
               ENDIF
               cos_small_angle = dcos(STL_SMALL_ANGLE / 180.0 * PI)

               IF(DABS(cos_angle)>cos_small_angle) THEN
                  IGNORE_CURRENT_FACET = .TRUE.    ! Ignore small facets
               ENDIF

               NORM = sqrt(N1**2+N2**2+N3**2)

               IF(NORM>TOL_STL) THEN
                  N1 = N1 /NORM
                  N2 = N2 /NORM
                  N3 = N3 /NORM
               ELSE
                  IGNORE_CURRENT_FACET = .TRUE.  ! Ignore facets with zero normal vector
               ENDIF

               IF(IGNORE_CURRENT_FACET) THEN
                  IGNORED_FACETS = IGNORED_FACETS + 1

               ELSE                              ! Save vertex coordinates for valid facets
                  N_FACETS_DES = N_FACETS_DES + 1
                  STL_GROUP(N_STL_GROUP)%END = N_FACETS_DES
                  NORM_FACE(1,N_FACETS_DES) = N1
                  NORM_FACE(2,N_FACETS_DES) = N2
                  NORM_FACE(3,N_FACETS_DES) = N3
                  VERTEX(1,1,N_FACETS_DES) = SCALE_STL*V1x + TX_STL
                  VERTEX(1,2,N_FACETS_DES) = SCALE_STL*V1y + TY_STL
                  VERTEX(1,3,N_FACETS_DES) = SCALE_STL*V1z + TZ_STL
                  VERTEX(2,1,N_FACETS_DES) = SCALE_STL*V2x + TX_STL
                  VERTEX(2,2,N_FACETS_DES) = SCALE_STL*V2y + TY_STL
                  VERTEX(2,3,N_FACETS_DES) = SCALE_STL*V2z + TZ_STL
                  VERTEX(3,1,N_FACETS_DES) = SCALE_STL*V3x + TX_STL
                  VERTEX(3,2,N_FACETS_DES) = SCALE_STL*V3y + TY_STL
                  VERTEX(3,3,N_FACETS_DES) = SCALE_STL*V3z + TZ_STL

               ENDIF

            ENDIF
         ENDDO

930      WRITE(ERR_MSG,100) 'Done reading file.'
         CALL LOG_STATUS()
         CLOSE(333)
      ENDDO


      XMIN_STL = MINVAL(VERTEX(:,1,1:N_FACETS))
      XMAX_STL = MAXVAL(VERTEX(:,1,1:N_FACETS))
      YMIN_STL = MINVAL(VERTEX(:,2,1:N_FACETS))
      YMAX_STL = MAXVAL(VERTEX(:,2,1:N_FACETS))
      ZMIN_STL = MINVAL(VERTEX(:,3,1:N_FACETS))
      ZMAX_STL = MAXVAL(VERTEX(:,3,1:N_FACETS))

      WRITE(ERR_MSG,2000)'STL file(s) successfully read.'
      CALL LOG_STATUS()
      WRITE(ERR_MSG,110)'Total number of facets read =',N_FACETS + IGNORED_FACETS
      CALL LOG_STATUS()
      WRITE(ERR_MSG,110)'Number of valid facets      =',N_FACETS
      CALL LOG_STATUS()
      WRITE(ERR_MSG,110)'Number of ignored facets    =',IGNORED_FACETS
      CALL LOG_STATUS()
      WRITE(ERR_MSG,100)'Geometry range from stl file:'
      CALL LOG_STATUS()
      IF(SCALE_STL/=ONE) THEN
         WRITE(ERR_MSG,5000)'Scaling factor:',SCALE_STL
      ENDIF
      ABSTRANS = dabs(TX_STL)+dabs(TY_STL)+dabs(TZ_STL)
      IF(ABSTRANS>TOL_STL) THEN
         WRITE(ERR_MSG,3000)'Translation vector (X,Y,Z)=',TX_STL,TY_STL,TZ_STL
         CALL LOG_STATUS()
      ENDIF
      WRITE(ERR_MSG,4000)'x-range = ', XMIN_STL,XMAX_STL
      CALL LOG_STATUS()
      WRITE(ERR_MSG,4000)'y-range = ', YMIN_STL,YMAX_STL
      CALL LOG_STATUS()
      WRITE(ERR_MSG,4000)'z-range = ', ZMIN_STL,ZMAX_STL
      CALL LOG_STATUS()
      WRITE(ERR_MSG,4000)''
      CALL LOG_STATUS()

      XMIN_STL = XMIN_STL - 10.0*TOL_STL
      XMAX_STL = XMAX_STL + 10.0*TOL_STL
      YMIN_STL = YMIN_STL - 10.0*TOL_STL
      YMAX_STL = YMAX_STL + 10.0*TOL_STL
      ZMIN_STL = ZMIN_STL - 10.0*TOL_STL
      ZMAX_STL = ZMAX_STL + 10.0*TOL_STL

!      IF(XMIN_STL<ZERO) XMIN_STL=ZERO
!      IF(XMAX_STL>XLENGTH) XMAX_STL=XLENGTH
!      IF(YMIN_STL<ZERO) YMIN_STL=ZERO
!      IF(YMAX_STL>YLENGTH) YMAX_STL=YLENGTH
!      IF(ZMIN_STL<ZERO) ZMIN_STL=ZERO
!      IF(ZMAX_STL>ZLENGTH) ZMAX_STL=ZLENGTH




      DO NG=1,N_STL_GROUP
         WRITE(*,*) 'STL GROUP: ', NG
         WRITE(*,*) 'TYPE     : ', STL_GROUP(NG)%TYPE
         WRITE(*,*) 'FILENAME : ', STL_GROUP(NG)%FILENAME
         WRITE(*,*) 'ID       : ', STL_GROUP(NG)%ID
         WRITE(*,*) 'START    : ', STL_GROUP(NG)%START
         WRITE(*,*) 'END      : ', STL_GROUP(NG)%END
      ENDDO


      RETURN
! HERE IF AN ERROR OCCURRED OPENING/READING THE FILE

 910  CONTINUE
      WRITE (ERR_MSG, 1500)
      call log_error()
 920  CONTINUE
      WRITE (ERR_MSG, 1600)
      call log_error()
! 930  CONTINUE
!     WRITE (ERR_MSG, 1700)
!     call log_error()

 1500 FORMAT(/1X,70('*')//' From: GET_STL_DATA_IS',/' Message: ',&
      'Unable to open stl file',/1X,70('*')/)
 1600 FORMAT(/1X,70('*')//' From: GET_STL_DATA_IS',/' Message: ',&
      'Error while reading stl file',/1X,70('*')/)
 2000 FORMAT(1X,A)
 3000 FORMAT(1X,A,'(',F10.4,';',F10.4,';',F10.4,')')
 4000 FORMAT(1X,A,F10.4,' to ',F10.4)
 5000 FORMAT(1X,A,F10.4)
      END SUBROUTINE GET_STL_DATA_IS


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: TRANSLATE_STL                                          C
!  Purpose: Translates a group of STL facets                           C
!                                                                      C
!  Author: Jeff Dietiker                              Date: 09-JUN-20  C
!  Reviewer:                                          Date: **-***-**  C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      subroutine Translate_STL(Group, Translation)

      INTEGER, intent(in) :: Group
      DOUBLE PRECISION, DIMENSION(3), intent(in) :: Translation
      INTEGER :: N1, N2
      INTEGER :: V, N

         N1 = STL_GROUP(Group)%START
         N2 = STL_GROUP(Group)%END

         DO N=N1,N2
            DO V=1,3
               VERTEX(V,:,N) = VERTEX_INIT(V,:,N) + Translation(:)
            ENDDO
         ENDDO

      RETURN
      End subroutine Translate_STL

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: ROTATE_STL                                             C
!  Purpose: Rotates a group of STL facets                              C
!                                                                      C
!  Author: Jeff Dietiker                              Date: 09-JUN-20  C
!  Reviewer:                                          Date: **-***-**  C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      subroutine Rotate_STL(Group, Axis, Center, Theta)

      INTEGER, intent(in) :: Group
      CHARACTER(LEN=1), intent(in) :: Axis
      DOUBLE PRECISION, DIMENSION(3), intent(in) :: Center
      DOUBLE PRECISION, intent(in) :: Theta ! Rotation angle (radians)
      INTEGER :: N1, N2
      DOUBLE PRECISION, DIMENSION(3,3) :: R
      DOUBLE PRECISION :: Theta_deg
      INTEGER :: V, N

! Converts angle from radians to degrees
      Theta_deg = Theta / acos(-1.0d0) * 180.0

! Build rotation matrix
      SELECT CASE(Axis)
         CASE('x','X')
            CALL BUILD_X_ROTATION_MATRIX(Theta_deg, R)

         CASE('y','Y')
            CALL BUILD_Y_ROTATION_MATRIX(Theta_deg, R)

         CASE('z','Z')
            CALL BUILD_Z_ROTATION_MATRIX(Theta_deg, R)

         END SELECT

! Apply rotation
         N1 = STL_GROUP(Group)%START
         N2 = STL_GROUP(Group)%END
         DO N=N1,N2
            DO V=1,3
               VERTEX(V,:,N) = MATMUL(R, VERTEX_INIT(V,:,N) - Center(:)) + Center(:)
            ENDDO
            NORM_FACE(:,N) = MATMUL(R, NORM_FACE_INIT(:,N))
         ENDDO

      RETURN
      End subroutine Rotate_STL

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: COPY_INIT_STL_DATA                                     C
!  Purpose: Keep a copy f the initial stl files. Operations (such as   C
!           translation and rotations are applied to the original stl  C
!                                                                      C
!  Author: Jeff Dietiker                              Date: 09-JUN-20  C
!  Reviewer:                                          Date: **-***-**  C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE COPY_INIT_STL_DATA
      IMPLICIT NONE
      INTEGER :: NN, NG, VV

      IF(.NOT.ALLOCATED(VERTEX_INIT)) ALLOCATE(VERTEX_INIT(3,3,DIM_STL))
      IF(.NOT.ALLOCATED(NORM_FACE_INIT)) ALLOCATE(NORM_FACE_INIT(3,DIM_STL))

      DO NG=1,N_STL_GROUP
         IF(STL_GROUP(NG)%TYPE=='IS_STL') THEN
            do NN=STL_GROUP(NG)%START,STL_GROUP(NG)%END
               do VV = 1,3
                  VERTEX_INIT(VV,1:3,NN) = VERTEX(VV,1:3,NN)
               enddo
            NORM_FACE_INIT(1:3,NN) = NORM_FACE(1:3,NN)
            enddo
         ENDIF
      ENDDO
      RETURN
      END SUBROUTINE COPY_INIT_STL_DATA

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: GET_STL_GROUP_ID                                       C
!  Purpose: Given an stl file name (say 'is_0002.stl', return the      C
!           corresponding STL group ID                                 C
!                                                                      C
!  Author: Jeff Dietiker                              Date: 09-JUN-20  C
!  Reviewer:                                          Date: **-***-**  C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      INTEGER FUNCTION GET_GROUP_ID(STL_FILE)
      IMPLICIT NONE
      INTEGER :: NN, NG, VV
      CHARACTER(LEN=*), INTENT(in) :: STL_FILE

      GET_GROUP_ID = 0
      DO NG=1,N_STL_GROUP
         IF(TRIM(STL_GROUP(NG)%FILENAME)==STL_FILE) GET_GROUP_ID = NG
      ENDDO
      ! print*,'Group ID of ', STL_FILE , GET_GROUP_ID
      IF(GET_GROUP_ID == 0) THEN
         WRITE(ERR_MSG,*)  'ERROR: STL GROUP ID NOT FOUND FOR: ',STL_FILE
         CALL LOG_ERROR()
      ENDIF
      RETURN
      END FUNCTION GET_GROUP_ID

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: ASSIGN_BC_ID_TO_STL_GROUP                              C
!  Purpose: Given an STL group ID, assign a BC_ID to the list of       C
!           facets in this group. Needed for IS STL files if we want   C
!           to define a wall velocity.                                 C
!                                                                      C
!  Author: Jeff Dietiker                              Date: 09-JUN-20  C
!  Reviewer:                                          Date: **-***-**  C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE ASSIGN_BC_ID_TO_STL_GROUP(GROUP_ID, BC_ID)
      IMPLICIT NONE
      INTEGER :: NN, GROUP_ID, BC_ID

      DO NN = STL_GROUP(GROUP_ID)%START, STL_GROUP(GROUP_ID)%END
         BC_ID_STL_FACE(NN) = BC_ID
      ENDDO
      RETURN
      END SUBROUTINE ASSIGN_BC_ID_TO_STL_GROUP

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: DETECT STL_ALONG_MFIX_BOX                              C
!  Purpose: detects stl facets that are along existing MFiX domain box.C
!           This will help get better intersection                     C
!                                                                      C
!  Author: Jeff Dietiker                              Date: 27-SEPT-19 C
!  Reviewer:                                          Date: **-***-**  C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE DETECT_STL_ALONG_MFIX_BOX(NUMBER_OF_CGBC_PATCHES)

      IMPLICIT NONE
      INTEGER, INTENT(IN) :: NUMBER_OF_CGBC_PATCHES
      INTEGER :: FACET, DIR, V, CDIR, F
      INTEGER :: BC_COUNTER, BC_HOLDER
      DOUBLE PRECISION :: TOL_C, TOL_N, CVAL, VMEAN
      character(len=1), dimension(6), parameter :: dir_to_char = (/'e','w','n','s','t','b'/)

      INTEGER :: NN,IGNORED_FACETS
      LOGICAL :: PRESENT,KEEP_READING,IGNORE_CURRENT_FACET
      DOUBLE PRECISION ::v1x,v1y,v1z
      DOUBLE PRECISION ::v2x,v2y,v2z
      DOUBLE PRECISION ::v3x,v3y,v3z
      DOUBLE PRECISION ::x12,y12,z12,x13,y13,z13,dp,d12,d13
      DOUBLE PRECISION ::cos_angle,cos_small_angle
      DOUBLE PRECISION ::n1,n2,n3,NORM
      DOUBLE PRECISION ::ABSTRANS
      CHARACTER(LEN=32) ::TEST_CHAR,BUFF_CHAR
      CHARACTER(LEN=255) ::geometryfile(0:DIMENSION_BC)
      INTEGER :: BCV,NUMBER_OF_GEOMETRY_FILES,NUMBER_OF_BC_PATCHES
      LOGICAL :: BC_PATCH_FOUND_IN_STL(DIMENSION_BC)
      CHARACTER(LEN=100) :: FNAME
      integer :: stl_unit, nf
      LOGICAL :: MULTISTL
!---------------------------------------------------------------------

      type box_boundary
         double precision  :: uv_to_fluid(3) ! Unit vector pointing towards fluid
         character(len=10) :: plane_loc      ! Short description of plane location
         integer           :: c_dir          ! Coordinate to check (1=x, 2=y, 3=z)
         double precision  :: c_value        ! coordinate value
      end type box_boundary

      type(box_boundary), dimension(6) :: box_plane


      type bc_structure
         integer              :: bc_id          ! Boundary condition ID
         integer              :: bc_type        ! Boundary condition Type
         integer              :: n_facets       ! Number of facets
         integer, dimension(6):: n_f_dir        ! Number of facets along each box plane (1=e, 2=w, 3=n, 4=s, 5=t, 6=b)
         integer, allocatable :: facet_list(:)  ! list of facets
         logical              :: belongs_to_box ! .True. if bc patch is aligned
                                                !  with one of the MFiX bounding box
         character(len=1)     :: face_char      ! face (e/w/n/s/t/b)
         integer              :: face           ! face (1/2/3/4/5/6)
         double precision     :: x_e, x_w, y_n, y_s, z_t, z_b

      end type bc_structure


      type(bc_structure), allocatable :: bc_patch(:)


! DETERMINE WHICH BOUNDARY CONDITIONS NEED STL FILE
      WRITE(ERR_MSG,100)'From Detect_stl_along_mfix_box:'// new_line('') // &
                        'Grouping facets by BC patch...'

      CALL LOG_STATUS()

! Tolerance
      TOL_C = 1.0D-6   ! used to compare coordinates (x,y and z)
      TOL_N = 1.0D-6   ! used to compare normal vectors


! Define MFiX bounding box

! East side
      box_plane(1)%uv_to_fluid = (/-ONE, ZERO, ZERO/)
      box_plane(1)%plane_loc   = 'East'
      box_plane(1)%c_dir       = 1
      box_plane(1)%c_value     = x_max

! West side
      box_plane(2)%uv_to_fluid = (/ ONE, ZERO, ZERO/)
      box_plane(2)%plane_loc   = 'West'
      box_plane(2)%c_dir       = 1
      box_plane(2)%c_value     = x_min

! North side
      box_plane(3)%uv_to_fluid = (/ZERO, -ONE, ZERO/)
      box_plane(3)%plane_loc   = 'North'
      box_plane(3)%c_dir       = 2
      box_plane(3)%c_value     = y_max

! South side
      box_plane(4)%uv_to_fluid = (/ZERO,  ONE, ZERO/)
      box_plane(4)%plane_loc   = 'South'
      box_plane(4)%c_dir       = 2
      box_plane(4)%c_value     = y_min

! Top side
      box_plane(5)%uv_to_fluid = (/ZERO, ZERO, -ONE/)
      box_plane(5)%plane_loc   = 'Top'
      box_plane(5)%c_dir       = 3
      box_plane(5)%c_value     = z_max

! Bottom side
      box_plane(6)%uv_to_fluid = (/ZERO, ZERO,  ONE/)
      box_plane(6)%plane_loc   = 'Bottom'
      box_plane(6)%c_dir       = 3
      box_plane(6)%c_value     = z_min


! Initialize Boundary condition patches

      if(.not.allocated(bc_patch)) allocate(bc_patch(NUMBER_OF_CGBC_PATCHES))

      do bcv = 1,NUMBER_OF_CGBC_PATCHES
         if(.not.allocated(bc_patch(bcv)%facet_list)) allocate (bc_patch(bcv)%facet_list(N_FACETS))
         bc_patch(bcv)%n_facets = 0
         bc_patch(bcv)%n_f_dir(:) = 0
         bc_patch(bcv)%belongs_to_box = .false.
         bc_patch(bcv)%face = 0
         bc_patch(bcv)%face_char = ' '
         bc_patch(bcv)%x_e = -UNDEFINED
         bc_patch(bcv)%x_w =  UNDEFINED
         bc_patch(bcv)%y_n = -UNDEFINED
         bc_patch(bcv)%y_s =  UNDEFINED
         bc_patch(bcv)%z_t = -UNDEFINED
         bc_patch(bcv)%z_b =  UNDEFINED
      enddo



      STL_FACET_ALONG_BOX(:) = .FALSE.

      BC_COUNTER = 0
      BC_HOLDER  = 0  ! used to increment BC_COUNTER

! STL facets where read sequentially and are already grouped by increasing
! BC ID. bc_patch stores the list of facets. This may be overkill since
! storing the starting and ending index for each BC could be sufficient.


      DO FACET = 1, N_FACETS
! Identify BC counter
         IF(BC_ID_STL_FACE(FACET)>BC_HOLDER) THEN
            BC_COUNTER = BC_COUNTER + 1
            BC_HOLDER  = BC_ID_STL_FACE(FACET)
         ENDIF

! Store list of facets
         bc_patch(bc_counter)%bc_id = bc_holder
         bc_patch(bc_counter)%n_facets = bc_patch(bc_counter)%n_facets + 1
         bc_patch(bc_counter)%facet_list(bc_patch(bc_counter)%n_facets) = FACET

! Snap verices along MFiX bounding box to eliminate lack of precision in the STL
! file
         DO DIR=1,6 ! Repeat for each side of the box
            CDIR = box_plane(DIR)%c_dir
            CVAL = box_plane(DIR)%c_value
            DO V = 1,3 ! Test the three vertices of each facet
               IF(DABS(VERTEX(V,CDIR,FACET)-CVAL)<TOL_C) VERTEX(V,CDIR,FACET) = CVAL
            ENDDO
! Check if the facet is parallel to one of the box side
            DP = DOT_PRODUCT(box_plane(DIR)%uv_to_fluid,NORM_FACE(:,FACET))
            IF(DP>ONE-TOL_N) THEN
! If it is parallel to a side, check if the facet lies along the edge
! Note: Checking one vertex instead of th mean may be sufficient
               VMEAN = ZERO
               DO V = 1,3
                  VMEAN = VMEAN + VERTEX(V,CDIR,FACET)
               ENDDO
               VMEAN = VMEAN / 3.0D0
! Add facet to the list of facets belonging to a side
               IF(DABS(VMEAN-CVAL)<TOL_C) THEN
                  STL_FACET_ALONG_BOX(FACET) = .TRUE.
                  bc_patch(bc_counter)%n_f_dir(dir) = bc_patch(bc_counter)%n_f_dir(dir) + 1
               ENDIF
            ENDIF
         ENDDO  ! Loop of MFiX Box faces (e/w/n/s/t/b)

      ENDDO ! Loop over facets


! An STL file is said to belong to a side if all the facets belongs to a
! side. Record which side of the box
      do bcv = 1,bc_counter
         do dir = 1,6
            if(bc_patch(bcv)%n_f_dir(dir)==bc_patch(bcv)%n_facets) then
               bc_patch(bcv)%belongs_to_box = .true.
               bc_patch(bcv)%face = dir
               bc_patch(bcv)%face_char = dir_to_char(dir)
            endif
         enddo
! Identify a rectandular region encompassing the STL file
! This should be a plane
         do f = 1,bc_patch(bcv)%n_facets
            FACET = bc_patch(bcv)%facet_list(f)
            DO V = 1,3
               if(VERTEX(V,1,FACET)>bc_patch(bcv)%x_e) bc_patch(bcv)%x_e = VERTEX(V,1,FACET)
               if(VERTEX(V,1,FACET)<bc_patch(bcv)%x_w) bc_patch(bcv)%x_w = VERTEX(V,1,FACET)
               if(VERTEX(V,2,FACET)>bc_patch(bcv)%y_n) bc_patch(bcv)%y_n = VERTEX(V,2,FACET)
               if(VERTEX(V,2,FACET)<bc_patch(bcv)%y_s) bc_patch(bcv)%y_s = VERTEX(V,2,FACET)
               if(VERTEX(V,3,FACET)>bc_patch(bcv)%z_t) bc_patch(bcv)%z_t = VERTEX(V,3,FACET)
               if(VERTEX(V,3,FACET)<bc_patch(bcv)%z_b) bc_patch(bcv)%z_b = VERTEX(V,3,FACET)
            ENDDO
         enddo
! Convert the 'CG_###' BC type into a regular '###' BC type
         if(bc_patch(bcv)%belongs_to_box) then
            WRITE(ERR_MSG,110)'CG BC: ',bc_patch(bcv)%bc_id, ' is flush with MFiX box.'// new_line('') // &
                              'It will be converted to regular BC.'// new_line('') // &
                              'Original BC type : ', bc_type(bc_patch(bcv)%bc_id),new_line('') // &
                              'New BC type      : ', bc_type(bc_patch(bcv)%bc_id)(4:6)
            CALL LOG_INFO()

            bc_type(bc_patch(bcv)%bc_id)=bc_type(bc_patch(bcv)%bc_id)(4:6)
            bc_x_w(bc_patch(bcv)%bc_id) = bc_patch(bcv)%x_w
            bc_x_e(bc_patch(bcv)%bc_id) = bc_patch(bcv)%x_e
            bc_y_s(bc_patch(bcv)%bc_id) = bc_patch(bcv)%y_s
            bc_y_n(bc_patch(bcv)%bc_id) = bc_patch(bcv)%y_n
            bc_z_b(bc_patch(bcv)%bc_id) = bc_patch(bcv)%z_b
            bc_z_t(bc_patch(bcv)%bc_id) = bc_patch(bcv)%z_t
         endif
      enddo

! Remove the facets that were identified as belonging to the MFiX bounding box.
! These will not be used, and a regular BC will be applied instead.
      N_FACETS = 0
      STL_FACET_ALONG_BOX(:) = .false.
      do bcv = 1,bc_counter
         if(bc_patch(bcv)%belongs_to_box) then ! Remove facets
            bc_patch(bcv)%n_facets = 0
            bc_patch(bcv)%facet_list(:) = -1
         else
! We only need to "slide" the facets indices
! No need to make a copy, the new index will alays be smaller or equal to the
! old one, so there is no risk of overwriting any data
            do f = 1,bc_patch(bcv)%n_facets
               FACET = bc_patch(bcv)%facet_list(f)
               N_FACETS = N_FACETS + 1
               bc_patch(bcv)%facet_list(f) = N_FACETS
!               STL_FACET_ALONG_BOX(N_FACETS) = .FALSE.
                  DO V = 1,3
                     VERTEX(V,1:3,N_FACETS) = VERTEX(V,1:3,FACET)
                  ENDDO
                  NORM_FACE(1:3,N_FACETS)  = NORM_FACE(1:3,FACET)
                  BC_ID_STL_FACE(N_FACETS) = BC_ID_STL_FACE(FACET)
            enddo
         endif
         if(.FALSE.) THEN  ! Print debug info
            print*,'BC counter       = ',bcv
            print*,'BC ID            = ',bc_patch(bcv)%bc_id
            print*,'BC TYPE          = ',BC_TYPE(bc_patch(bcv)%bc_id)
            print*,'total # facets   = ',bc_patch(bcv)%n_facets
            print*,'faces list       = ',bc_patch(bcv)%facet_list(1:bc_patch(bcv)%n_facets)
            print*,'# faces per side = ',bc_patch(bcv)%n_f_dir(:)
            print*,' Belongs to box  = ',bc_patch(bcv)%belongs_to_box
            print*,' Box face        = ',bc_patch(bcv)%face
            print*,' Box face char   = ',bc_patch(bcv)%face_char
            print*,' BC extents:       '
            print*,' x-west          = ',bc_patch(bcv)%x_w
            print*,' x-east          = ',bc_patch(bcv)%x_e
            print*,' y-south         = ',bc_patch(bcv)%y_s
            print*,' y-north         = ',bc_patch(bcv)%y_n
            print*,' z-bottom        = ',bc_patch(bcv)%z_b
            print*,' z-top           = ',bc_patch(bcv)%z_t
         endif
      enddo

      WRITE(ERR_MSG,120)'Number of active facets (BCs not flush with MFiX box) = ', N_FACETS
      CALL LOG_INFO()

100  FORMAT(1X,A)
110  FORMAT(1X,A,I4,A,A,A,A)
120  FORMAT(1X,A,I6)


      END SUBROUTINE DETECT_STL_ALONG_MFIX_BOX


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: EVAL_STL_FCT_AT                                        C
!  Purpose: Assigns a value to f_stl at any NODE of cell IJK           C
!                                                                      C
!  Author: Jeff Dietiker                              Date: 30-JAN-09  C
!  Reviewer:                                          Date: **-***-**  C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE EVAL_STL_FCT_AT(TYPE_OF_CELL,IJK,NODE,f_stl,CLIP_FLAG,BCID)

      IMPLICIT NONE
      INTEGER :: BCID
      LOGICAL :: CLIP_FLAG
      DOUBLE PRECISION :: f_stl
      CHARACTER (LEN=*) :: TYPE_OF_CELL
      INTEGER :: IJK,IJKC,NODE
!---------------------------------------------------------------------

      IF(NODE==15) print*,'Warning: eval stl at node 15'

      IF(N_FACETS < 1) RETURN

      f_stl = UNDEFINED
      IJKC = IJK_OF_NODE(NODE)
      IF(IJKC<1.OR.IJKC>DIMENSION_3) THEN

!         print*,'myPE        =',myPE
!         print*,'DIMENSION_3 =',DIMENSION_3
!         print*,'IJKC        =',IJKC
!         print*,'IJK         =',IJK
!         print*,'I,J,K       =',I_OF(IJK),J_OF(IJK),K_OF(IJK)
!         print*,'NODE        =',NODE

!     Leave f_stl as undefined if we are out of bounds
!     This can happen at I=1 and NODE=4 for example
         RETURN
      ENDIF

!!      IF(NODE/=15.AND.NODE/=0) THEN
!!         IF(SNAP(IJKC)) THEN
!!            f_stl = ZERO
!!!            print*,'stlfc snapped=',IJKC,NODE,F_AT(IJKC)
!!            RETURN
!!         ENDIF
!!      ENDIF

      f_stl = F_AT(IJKC)
      RETURN

      CLIP_FLAG = .TRUE.
      RETURN

      END SUBROUTINE EVAL_STL_FCT_AT


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: intersect_line_with_facet                              C
!  Purpose: Finds the intersection between a facet                     C
!           and the line (xa,ya,za) and (xb,yb,zb).                    C
!                                                                      C
!  Author: Jeff Dietiker                              Date: 21-Feb-08  C
!  Reviewer:                                          Date:            C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

! Notes:
! This subroutine tries to find the intersection of a line AB with one
! of the facets defined in the stl file
! 1) Verify that intersection is possible. Only facets and line AB that
!    are not parallel are acceptable candidates
! 2) Find intersection point P of line AB with plane containing the facet.
!    Valid intersection point must be between A and B to continue.
! 3) Verify that point P is inside triangular facet

  SUBROUTINE INTERSECT_LINE_WITH_FACET(xa,ya,za,xb,yb,zb,FACET,INTERSECT_FLAG,xc,yc,zc)

      IMPLICIT NONE
      DOUBLE PRECISION:: xa,ya,za,xb,yb,zb,xc,yc,zc
      INTEGER :: FACET
      DOUBLE PRECISION:: NFx,NFy,NFz,nabx,naby,nabz
      DOUBLE PRECISION :: dot_denom,dot_num
      DOUBLE PRECISION :: VP1Ax,VP1Ay,VP1Az
      DOUBLE PRECISION :: Px,Py,Pz
      DOUBLE PRECISION :: tt
      DOUBLE PRECISION :: d_ac,d_bc
      LOGICAL :: INSIDE_FACET,INTERSECT_FLAG
!---------------------------------------------------------------------

      INTERSECT_FLAG = .FALSE.

! Facet normal vector (normalized)
      NFx = NORM_FACE(1,FACET)
      NFy = NORM_FACE(2,FACET)
      NFz = NORM_FACE(3,FACET)

! AB vector (NOT normalized)
      nabx = xb - xa
      naby = yb - ya
      nabz = zb - za

      dot_denom = NFx*nabx + NFy*naby + NFz*nabz

! 1) Verify that intersection is possible. Only facets and line AB that are
!    not parallel are acceptable candidates
      IF(dabs(dot_denom)<TOL_STL) THEN
         INTERSECT_FLAG = .FALSE.              ! No intersection (facet nornal
                                               ! and AB are perpendicular
         RETURN
      ELSE

! 2) Find intersection point P of line AB with plane containing the facet
!    Line AB is parametrized with parameter t (from t=0 at A to t=1 at B)
         VP1Ax = VERTEX(1,1,FACET) - xa
         VP1Ay = VERTEX(1,2,FACET) - ya
         VP1Az = VERTEX(1,3,FACET) - za
         dot_num = NFx*VP1Ax + NFy*VP1Ay + NFz*VP1Az
         tt = dot_num / dot_denom

! 3) Verify that point P is inside triangular facet
!    Facet is parametrized with (u,v), a point is inside triangle if:
!    u   >= 0  , and
!    v   >= 0  , and
!    u+v <= 1

         IF((tt>=ZERO).AND.(tt<=ONE)) THEN     ! Intersection between A and B
! Now test if intersection point is inside triangle
            Px = xa + tt*nabx
            Py = ya + tt*naby
            Pz = za + tt*nabz
            CALL IS_POINT_INSIDE_FACET(Px,Py,Pz,FACET,INSIDE_FACET)
            IF(INSIDE_FACET) THEN
               INTERSECT_FLAG = .TRUE.              ! Valid intersection
               xc = Px                              ! point P is inside triangle)
               yc = Py
               zc = Pz
            ELSE
               INTERSECT_FLAG = .FALSE.             ! Invalid intersection
               RETURN                               ! point P is outside triangle)
            ENDIF
         ELSE
            INTERSECT_FLAG = .FALSE.           ! No intersection between A and B
            RETURN
         ENDIF
      ENDIF

      d_ac = sqrt((xc-xa)**2 + (yc-ya)**2 + (zc-za)**2)
      d_bc = sqrt((xc-xb)**2 + (yc-yb)**2 + (zc-zb)**2)

      IF(d_ac<TOL_STL.OR.d_bc<TOL_STL) THEN
         INTERSECT_FLAG = .FALSE.             ! Exclude corner intersection
      ENDIF

      RETURN
      END SUBROUTINE INTERSECT_LINE_WITH_FACET


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: IS_POINT_INSIDE_FACET                                  C
!  Purpose: Verifies that point P is inside facet                      C
!                                                                      C
!  Author: Jeff Dietiker                              Date: 21-Feb-08  C
!  Reviewer:                                          Date:            C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
  SUBROUTINE IS_POINT_INSIDE_FACET(Px,Py,Pz,FACET,INSIDE_FACET)

      IMPLICIT NONE

      INTEGER :: FACET,VV
      DOUBLE PRECISION :: NFx,NFy,NFz
      DOUBLE PRECISION :: Px,Py,Pz,V0x,V0y,V0z,V1x,V1y,V1z,V2x,V2y,V2z
      DOUBLE PRECISION :: V2n,V2nx,V2ny,V2nz
      DOUBLE PRECISION :: dot00,dot01,dot02,dot11,dot12,dot_check
      DOUBLE PRECISION :: Inv_denom
      DOUBLE PRECISION :: u,v
      LOGICAL :: U_POSITIVE,V_POSITIVE,UPVL1,INSIDE_FACET

      DOUBLE PRECISION :: Vx,Vy,Vz
      DOUBLE PRECISION :: D(3),MINVAL_D,DH
!---------------------------------------------------------------------

! If point P is very close to one of the Facet vertices, consider that
! P belongs to facet, and return.
         DO VV = 1,3
            Vx = VERTEX(VV,1,FACET)
            Vy = VERTEX(VV,2,FACET)
            Vz = VERTEX(VV,3,FACET)
            D(VV) = sqrt((Px - Vx)**2 + (Py - Vy)**2 + (Pz - Vz)**2 )
         ENDDO
         MINVAL_D = MINVAL(D)
         IF(MINVAL_D < TOL_STL) THEN
            INSIDE_FACET = .TRUE.
            RETURN
         ENDIF


! Facet normal vector (normalized)
      NFx = NORM_FACE(1,FACET)
      NFy = NORM_FACE(2,FACET)
      NFz = NORM_FACE(3,FACET)

! This subroutine verifies that point P is inside triangular facet
! Facet is parametrized with (u,v), a point is inside triangle if:
! u   >= 0  , and
! v   >= 0  , and
! u+v <= 1
      V0x = VERTEX(2,1,FACET) - VERTEX(1,1,FACET)
      V0y = VERTEX(2,2,FACET) - VERTEX(1,2,FACET)
      V0z = VERTEX(2,3,FACET) - VERTEX(1,3,FACET)
      V1x = VERTEX(3,1,FACET) - VERTEX(1,1,FACET)
      V1y = VERTEX(3,2,FACET) - VERTEX(1,2,FACET)
      V1z = VERTEX(3,3,FACET) - VERTEX(1,3,FACET)
      V2x = Px - VERTEX(1,1,FACET)
      V2y = Py - VERTEX(1,2,FACET)
      V2z = Pz - VERTEX(1,3,FACET)
      dot_check = NFx*V2x + NFy*V2y + NFz*V2z

! Normalize V2 before dot product
      V2n = sqrt(V2x**2 + V2y**2 + V2z**2 )
      IF(V2n>ZERO) THEN
         V2nx = V2x/V2n
         V2ny = V2y/V2n
         V2nz = V2z/V2n
      ELSE
         V2nx = V2x
         V2ny = V2y
         V2nz = V2z
      ENDIF
      DH = NFx * V2nx + NFy * V2ny + NFz * V2nz
!      DH = NFx * V2x + NFy * V2y + NFz * V2z

      IF(dabs(DH)>TOL_STL_DP) THEN          ! reject points that do not
         INSIDE_FACET = .FALSE.             ! belong to plane containing facet
         RETURN
      ENDIF

      dot00 = V0x*V0x + V0y*V0y + V0z*V0z
      dot01 = V0x*V1x + V0y*V1y + V0z*V1z
      dot02 = V0x*V2x + V0y*V2y + V0z*V2z
      dot11 = V1x*V1x + V1y*V1y + V1z*V1z
      dot12 = V1x*V2x + V1y*V2y + V1z*V2z
      Inv_denom = ONE / (dot00*dot11 - dot01*dot01)
      u = (dot11*dot02 - dot01*dot12) * Inv_denom
      v = (dot00*dot12 - dot01*dot02) * Inv_denom
      U_POSITIVE = (u>=-TOL_STL_DP)
      V_POSITIVE = (v>=-TOL_STL_DP)
      UPVL1 = ((u+v)<=ONE+TOL_STL_DP)
      INSIDE_FACET = (U_POSITIVE.AND.V_POSITIVE.AND.UPVL1)

      RETURN

   END SUBROUTINE IS_POINT_INSIDE_FACET

END MODULE GET_STL_DATA_MOD
